__host__ __device__ float EvalGenAlgSurf(float x, float y, float z)
{
	//return x*x+y*y+z*z-1.0f;										//0 sphere
	//return 4*(x*x+y*y+z*z)-16*x*y*z;								//1 cayley
	//return -5*(x*x*(y+z)+y*y*(x+z)+z*z*(y+x))+2*(x*y+x*z+y*z);		//2 cayley
	return x*x+y*y-z*z*z*z*(1-z);									//3 kiss
	//return (1+x*x)*(1+y*y)*(1+z*z)+8*x*y*z-2.0;					//4 McMullen K3 model
	//return (x*x+3)*(y*y+3)*(z*z+3)-32*(x*y*z+1);					//5 Kusner-Schmitt
	//float x2=x*x,y2=y*y,z2=z*z; return x2*(x2-5)+y2*(y2-5)+z2*(z2-5)+11.8;				//6 Tangle
	//float t=(x*x+9.0f/4*y*y+z*z-1);return t*t*t-z*z*z*(x*x+9.0f/80*y*y);//7 Heart
	//float x2=x*x,y2=y*y,z2=z*z,t=(x2+y2+z2-2*y-1);return (t+4*y)*(t*t-8*z2)+16*x*z*t;	//8 Klein Bottle

	//float x2=x*x,y2=y*y,z2=z*z;return z2*(x2-z2)-2*z*(x2-z2)+x2-z2-(x2-z)*(x2-z)-y2*(y2+2*x2+z2-2*z-1);			//9 Nose 2
	//float x2=x*x,y2=y*y,z2=z*z;return z2*(x2-z2)-2*z*(x2-z)+x2-z2-(x2-z)*(x2-z)-y2*(y2+2*x2+z2-2*z-1);		//9 Nose

	//return 0.04*(y-1)+(x*x+z*z)*y;									//10 witch hat
	//float x2=x*x,y2=y*y,z2=z*z;return 4*x2*(x2+y2+z2+z)+y2*(y2+z2-1);	//11 Cross cap
	//float x2=x*x,y2=y*y,z2=z*z;return x2*y2+x2*z2+z2*y2+x*y*z;			//12 Steiner
	//float x2=x*x,y2=y*y,z2=z*z,t1=x2+y2+z2-13,t2=3*x2+y2-4*z2-12; return 4*t1*t1*t1+27*t2*t2;	//13 Hunt
	//float x2=x*x,y2=y*y,z2=z*z;return 4*x2*(x2+y2+z2)-y2*(1-y2-z2);			//14 Mitre
	//return y*y*y+z*z*z-6*y*z;													//15 Nodal cubic
	//float x2=x*x,y2=y*y,z2=z*z,t=2*x2+y2+z2-1; return t*t*t-x2*z2*z/10-y2*z2*z;	//16 Heart 2
	//return x*x*(x*x-1)+y*y+z*z;														//17 Lemniscate of Gerono
	//float x2=x*x,y2=y*y,z2=z*z,m=x2+y2-1,n=y2+z2-1,q=z2+x2-1,c1=0.075,c2=3;return (m*m+z2)*(n*n+x2)*(q*q+y2)-c1*c1*(1+c2*(x2+y2+z2)); //18 Orthocircle

	//float dis=(.3-sqrtf(x*x+y*y)); return dis*dis+z*z-.04;				//torus 
	//float x2=x*x,y2=y*y,z2=z*z,dis=x2+y2+z2+.09-.04; return dis*dis-4*.09*(x2+y2);//torus 4degree  19

	//float t=(1+sqrt(5.0))*.5,t2=(3+  sqrt(5.0))*.5,x2=x*x,y2=y*y,z2=z*z,w2=1; return 4*(t2*x2-y2)*(t2*y2-z2)*(t2*z2-x2)-(1+2*t)*(x2+y2+z2-w2)*(x2+y2+z2-w2)*w2; //20 barth sextic
	//float t=(1+sqrt(5.0))*.5,t4=(7+3*sqrt(5.0))*.5,x2=x*x,y2=y*y,z2=z*z,w2=1,s=(x2+y2+z2-w2)*(x2+y2+z2-(2-t)*w2); return 8*(x2-t4*y2)*(y2-t4*z2)*(z2-t4*x2)*(x2*x2+y2*y2+z2*z2-2*x2*y2-2*x2*z2-2*y2*z2)+(3+5*t)*s*s*w2;  //21 barth decic
	//float x2=x*x,y2=y*y,z2=z*z,yz2=y2+z2; return 4*(x2*x2+yz2*yz2)+17*x2*yz2-20*(x2+yz2)+17;  //22 torus Gumdrop
	//float x2=x*x,y2=y*y,z2=z*z; return 0.04-x2*x2*(1-2*x2+x2*x2)+y2*(2*x2*(1-x2)-y2)-z2;  //23 Double torus
	//float x2=x*x,y2=y*y,z2=z*z,c2=2.5,t=(x2+y2+1-c2); t=(t-x)*(t+x); t=t*t; return t+z2-.125;	//24 pretzel
	//return 8*z*z+(2*x+1)*(3*y*y-(x-1)*(x-1));					//25 Tritrumpet
	//float x2=x*x,t=(x2*(1-x2)-y*y); return  t*t+z*z-0.01;	//26 Bretzel
	//return x*x+y*y-z*z*(1-z);  //27 dingdong surface

	//return -(cos(x) + cos(y) + cos(z));  //28 k3DSurf schwartz
	//return (x*x+y*y-1)*(x*x+z*z-1)*(y*y+z*z-1)-1;		//29 TriCylinder
	//return   -(x*x+y*y+z*z)+cos(5*x)+cos(5*y)+cos(5*z)-.1;	//30 blobSchwartz
}

__host__ __device__ 
void EvalGenAlgSurfDf(float x, float y, float z, float df[3])
{
	//df[0]=2*x; df[1]=2*y; df[2]=2*z;										//0 sphere
	df[0]=2*x;df[1]=2*y;df[2]=z*z*z*(5*z-4);								//3 kiss
	//float x2=x*x,y2=y*y,z2=z*z; df[0]=(1+y2)*(1+z2)*2*x+8*y*z; df[1]=(1+x2)*(1+z2)*2*y+8*x*z;df[2]=(1+x2)*(1+y2)*2*z+8*x*y;			//4 McMullen K3 model
	//df[0]=x*(4*x*x-10); df[1]=y*(4*y*y-10); df[2]=z*(4*z*z-10);				//6 Tangle
	//float x2=x*x,y2=y*y,z2=z*z,t=(x2+y2+z2-2*y-1);df[0]=2*x*(t*t-8*z2)+(t+4*y)*(4*t*x)+16*z*(t+x2*2);df[1]=2*(y+1)*(t*t-8*z2)+(t+4*y)*(4*t*(y-1))+32*x*z*(y-1);df[2]=2*z*(t*t-8*z2)+(t+4*y)*(4*t*z-16*z)+16*x*(t+2*z2);	//8 Klein Bottle

	//df[0]=x*(4*x*x-2);df[1]=2*y; df[2]=2*z;										//17 Lemniscate of Gerono
	//float x2=x*x,y2=y*y; df[0]=4*x*((1-2*x2)*y2-x2*(1-3*x2+2*x2*x2)); df[1]=2*y*(3*x2*(1-x2)-2*y2); df[2]=-2*z; //23 Double torus
	//df[0]=2*x; df[1]=2*y; df[2]= 3*z*z-2*z;  //27 dingdong surface
	//df[0]= sin(x); df[1]=sin(y); df[2]=sin(z);  //28 k3DSurf schwartz
	//float x2=x*x,y2=y*y,z2=z*z,t=(x2+y2+z2-1),t0=t-z2,t1=t-y2,t2=t-x2; df[0]=2*x*(t0+t1)*t2,df[1]=2*y*(t0+t2)*t1,df[2]=2*z*(t1+t2)*t0;		//29 TriCylinder
	//df[0]=-2*x-5*sin(5*x),df[1]=-2*y-5*sin(5*y),df[2]=-2*z-5*sin(5*z);	//30 blobSchwartz


	//return x*x+y*y+z*z-1.0f;										//0 sphere
	//return 4*(x*x+y*y+z*z)-16*x*y*z;								//1 cayley
	//return -5*(x*x*(y+z)+y*y*(x+z)+z*z*(y+x))+2*(x*y+x*z+y*z);		//2 cayley
	//return (1+x*x)*(1+y*y)*(1+z*z)+8*x*y*z-2.1;					//4 McMullen K3 model
	//return (x*x+3)*(y*y+3)*(z*z+3)-32*(x*y*z+1);					//5 Kusner-Schmitt
	//float t=(x*x+9.0f/4*y*y+z*z-1);return t*t*t-z*z*z*(x*x+9.0f/80*y*y);//7 Heart

	//float x2=x*x,y2=y*y,z2=z*z;return z2*(x2-z2)-2*z*(x2-z2)+x2-z2-(x2-z)*(x2-z)-y2*(y2+2*x2+z2-2*z-1);			//9 Nose 2
	//float x2=x*x,y2=y*y,z2=z*z;return z2*(x2-z2)-2*z*(x2-z)+x2-z2-(x2-z)*(x2-z)-y2*(y2+2*x2+z2-2*z-1);		//9 Nose

	//return 0.04*(y-1)+(x*x+z*z)*y;									//10 witch hat
	//float x2=x*x,y2=y*y,z2=z*z;return 4*x2*(x2+y2+z2+z)+y2*(y2+z2-1);	//11 Cross cap
	//float x2=x*x,y2=y*y,z2=z*z;return x2*y2+x2*z2+z2*y2+x*y*z;			//12 Steiner
	//float x2=x*x,y2=y*y,z2=z*z,t1=x2+y2+z2-13,t2=3*x2+y2-4*z2-12; return 4*t1*t1*t1+27*t2*t2;	//13 Hunt
	//float x2=x*x,y2=y*y,z2=z*z;return 4*x2*(x2+y2+z2)-y2*(1-y2-z2);			//14 Mitre
	//return y*y*y+z*z*z-6*y*z;													//15 Nodal cubic
	//float x2=x*x,y2=y*y,z2=z*z,t=2*x2+y2+z2-1; return t*t*t-x2*z2*z/10-y2*z2*z;	//16 Heart 2
	//float x2=x*x,y2=y*y,z2=z*z,m=x2+y2-1,n=y2+z2-1,q=z2+x2-1,c1=0.075,c2=3;return (m*m+z2)*(n*n+x2)*(q*q+y2)-c1*c1*(1+c2*(x2+y2+z2)); //18 Orthocircle

	//float dis=(.3-sqrtf(x*x+y*y)); return dis*dis+z*z-.04;				//torus 
	//float x2=x*x,y2=y*y,z2=z*z,dis=x2+y2+z2+.09-.04; return dis*dis-4*.09*(x2+y2);//torus 4degree  19

	//float t=(1+sqrt(5.0))*.5,t2=(3+  sqrt(5.0))*.5,x2=x*x,y2=y*y,z2=z*z,w2=1; return 4*(t2*x2-y2)*(t2*y2-z2)*(t2*z2-x2)-(1+2*t)*(x2+y2+z2-w2)*(x2+y2+z2-w2)*w2; //20 barth sextic
	//float t=(1+sqrt(5.0))*.5,t4=(7+3*sqrt(5.0))*.5,x2=x*x,y2=y*y,z2=z*z,w2=1,s=(x2+y2+z2-w2)*(x2+y2+z2-(2-t)*w2); return 8*(x2-t4*y2)*(y2-t4*z2)*(z2-t4*x2)*(x2*x2+y2*y2+z2*z2-2*x2*y2-2*x2*z2-2*y2*z2)+(3+5*t)*s*s*w2;  //21 barth decic
	//float x2=x*x,y2=y*y,z2=z*z,yz2=y2+z2; return 4*(x2*x2+yz2*yz2)+17*x2*yz2-20*(x2+yz2)+17;  //22 torus Gumdrop

	//float x2=x*x,y2=y*y,z2=z*z,c2=2.5,t=(x2+y2+1-c2); t=(t-x)*(t+x); t=t*t; return t+z2-.125;	//24 pretzel
	//return 8*z*z+(2*x+1)*(3*y*y-(x-1)*(x-1));					//25 Tritrumpet
	//float x2=x*x,t=(x2*(1-x2)-y*y); return  t*t+z*z-0.01;	//26 Bretzel


}

__host__ __device__
float EvalGenAlgSurfPolarEye(float x, float y, float z, float e[3])
{	
	//return e[0]*2*x+e[1]*2*y+e[2]*2*z-2.0f;																							//0 sphere
	//return 4*(e[0]*(2*x-4*y*z)+e[1]*(2*y-4*x*z)+e[2]*(2*z-4*x*y)+x*x+y*y+z*z);															//1 cayley
	//float a=z+y,b=x+z,c=x+y,x2=x*x,y2=y*y,z2=z*z; return e[0]*(2*a-5*(y2+z2)-10*x*a)+e[1]*(2*b-5*(x2+z2)-10*y*b)+e[2]*(2*c-5*(x2+y2)-10*z*c)+2*(x*y+x*z+y*z);	//2 cayley
	//return e[0]*(2*x)+e[1]*(2*y)+e[2]*(z*z*z*(5*z-4))+3*(x*x+y*y)-z*z*z*z;									//3 kiss
	//float a=x*x+1,b=y*y+1,c=z*z+1;return e[0]*(2*x*b*c+8*y*z)+e[1]*(2*y*a*c+8*x*z)+e[2]*(2*z*a*b+8*x*y)+2*(b*c+a*c+a*b)+24*x*y*z-12.6;			//4 McMullen K3 model
	//float a=x*x+3,b=y*y+3,c=z*z+3; return e[0]*(2*x*b*c-32*y*z)+e[1]*(2*y*a*c-32*x*z)+e[2]*(2*z*a*b-32*x*y)+2*(b*c+a*c+a*b)-32*(3*x*y*z+6);					//5 Kusner-Schmitt
	float x2=x*x,y2=y*y,z2=z*z; return e[0]*(x*(4*x2-10))+e[1]*(y*(4*y2-10))+e[2]*(z*(4*z2-10))+4*11.8-10*(x2+y2+z2);			//6 Tangle
	//float x2=x*x,y2=y*y,z2=z*z,yy=2*y,t=(x2+y2+z2-2*y-1),r=t*t-8*z2; 	return e[0]*(2*x*r+4*x*t*(t+4*y)+16*z*(t+2*x2))+e[1]*((yy+2)*r+2*t*(t+4*y)*(yy-2)+16*x*z*(yy-2))+e[2]*(2*z*r+4*z*(t-4)*(t+4*y)+16*x*(t+2*z2))+(yy-2)*r+(t+4*y)*(2*t*(-yy-2)-16*z2)+16*x*z*(2*t-yy-2) ;	//8 Klein Bottle
	//return e[1]*(3*y*y-6*z)+e[2]*(3*z*z-6*y)-6*y*z;													//15 Nodal cubic
	//return e[0]*(4*x*x*x-2*x)+e[1]*(2*y)+e[2]*(2*z)+2*(y*y+z*z-x*x);									//17 Lemniscate of Gerono
	//float x2=x*x,y2=y*y,z2=z*z,dis=x2+y2+z2+.09-.04; return e[0]*(4*x*(dis-2*.09))+e[1]*(4*y*(dis-2*.09))+e[2]*(4*z*dis)+(4*(.09-.04)*dis-2*4*.09*(x2+y2));	//19 torus 4degree
	//return 1;//the f_w is hard //23 Double torus
	//float df0=2*x, df1=2*y, df2=3*z*z-2*z; return e[0]*df0+e[1]*df1+e[2]*df2+x*x+y*y-z*z;//27 dingdong surface

	//The rest are not computed.
}

__host__ __device__
float EvalGenAlgSurfPolar(float x, float y, float z)
{
	//float e[3]={4.0f, -6.2f, -4.0f};
	float e[3]={0.273231f, 2.519702f, 1.403071f};
	return EvalGenAlgSurfPolarEye(x, y, z, e);
}


